Spin-polarized low-density neutron matter 
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Low-density neutron matter is relevant to the study of neutron-rich nuclei and neutron star 
crusts. Unpolarized neutron matter has been studied extensively over a number of decades, while 
experimental guidance has recently started to emerge from the field of ultracold atomic gases. In 
this work, we study population-imbalanced neutron matter (possibly relevant to magnetars and to 
density functionals of nuclei) applying a Quantum Monte Carlo method that has proven successful 
in the field of cold atoms. We report on the first ab initio simulations of superfluid low-density 
polarized neutron matter. For systems with small imbalances, we find a linear dependence of the 
energy on the polarization, the proportionality coefficient being dependent on the density. We also 
present results for the momentum and pair distributions of the two fermionic components. 

PACS numbers: 21.65.-f, 03.75.Ss, 05.30.Fk, 26.60.-c 



I. INTRODUCTION 

The inner crust of a neutron star is widely considered 
to be composed of a lattice of neutron-rich nuclei along 
with a gas of neutrons and electrons. The gas of neutrons 
is expected to be superfluid at weak to intermediate cou- 
pling. Thus, low-density neutron matter is intrinsically 
connected to the strongly coupled fermion many-body 
problem, necessitating accurate calculations (or simula- 
tions). Low-density neutron matter is of relevance to the 
static and dynamic properties of the neutron star crust, 
which can lead to observable behavior. [Tj-Q Outside 
the observational realm, neutron matter computations 
also hold significance in the context of traditional nu- 
clear physics: equation of state results at densities close 
to the nuclear saturation density have been used for some 
time to constrain Skyrme and other density functional 
approaches to heavy nuclei, while the density-dependence 
of the 1 So gap in low-density neutron matter has also 
been used to constrain Skyrmc-Hartrec-Fock-Bogoliubov 
treatments in their description of neutron- rich nuclei Q ■ 
The potential significance of such calculations has led to 
a series of publications on the equation of state of low- 
density neutron matter over the last few decades. 

Parallel developments in a separate field of physics 
have recently provided new insight as well as the promise 
of direct experimental constraints: experiments with ul- 
tracold atomic gases of fermions are now being carried 
out in a number of labs around the world. In some cases, 
atomic gases in traps contain sufficiently many particles 
that the local-density approximation is valid. As a re- 
sult, such experiments can measure the energy fl7L [Tc| 
and pairing gap H(| of homogeneous strongly inter- 
acting matter. For cold fermionic atoms, the two-particle 
interaction can be directly tuned using a magnetic field 
through so-called Feshbach resonances to produce a spe- 
cific scattering length a, while the effective range r e be- 
tween the atoms is considerably smaller than the average 
interparticle distance, and thus essentially zero. These 
conditions are analogous to low-density neutron matter, 
where the particle-particle interaction has a scattering 



length which is considerable, ~ —18.5 fm, and is there- 
fore larger than the average interneutron spacing. On 
the other hand, cold atoms and low-density neutron mat- 
ter are clearly distinct systems: first, for neutron matter 
the effective range is much smaller than the scattering 
length, r e ks 2.7 fm, so \r e /a\ ~ 0.15, but only at very 
low densities is the effective range much smaller than the 
interparticle spacing. Second, the neutron-neutron (NN) 
interaction is not strictly limited to s- waves, implying a 
complicated spin dependence. Third, three-neutron in- 
teractions (NNN) are in principle also present. The last 
two points can be remedied if one studies relatively low- 
densities, i.e. at most an order of magnitude smaller than 
the nuclear saturation density. The first point might be 
addressed in the framework of cold atomic experiments 
in the future: it may be possible to use narrow and wide 
resonances in cold atoms to study this experimentally. [2l| 

The above discussion is limited to unpolarized neutron 
matter, i.e. to the case of two species of particles, conven- 
tionally called spin- up (t) and spin-down (4-), with equal 
populations. However, the general case of imbalanccd 
systems has also been studied extensively: limiting our- 
selves for the moment to neutron matter, calculations of 
the magnetic susceptibility have bee n ap pearing consis- 
tently since the discovery of pulsars. |22| Given that the 
magnitude of the pairing gap is approximately 1 MeV, a 
magnetic field in a neutron star crust would have to be 
larger than 10 16 — 10 17 G to polarize neutron matter. The 
question of spin-polarized neutron matter is thus in prin- 
ciple relevant to the objects known as magnetars, which 
have surface magnetic fields of 10 14 — 10 15 G. Thus, the 
promise of observational insight into these objects has led 
a number of theoretical groups to study spin-polarized 
neutron matter and the associated question of a possible 
ferromagnetic instability at large density. |23l - [29| 

In this work, we take a step back and address neu- 
tron matter with a finite spin-polarization (population 
imbalance) at low density. This system is in princi- 
ple relevant to neutron star observations: in a realistic 
neutron-star crust polarization may appear at lower den- 
sities than for infinite matter. For example, the spin- 
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orbit splitting around a large nucleus might help favor 
polarization at a lower magnetic field than would be re- 
quired for the bulk. Furthermore, our results could also 
be used as an input to or benchmark for phenomenolog- 
ical theories of terrestrial nuclei. On a different note, 
as already mentioned, the region of low-density is close 
to the physics of cold atoms, since there the NN inter- 
action is simpler and the NNN interaction is minimal. 
Thus, we can use approaches already verified in the lab- 
oratory with ultracold atomic gases. Importantly, we 
use a Quantum Monte Carlo (QMC) method which has 
been applied, in p revious works, to polarized cold atomic 
systems, (3(| HU as well as to unpolarized cold atoms 
and neutron matter [321433 ] . Our earlier works were con- 
sistent both with experimental measurements and with 
the analytically known behavior of the energy and the 
gap at vanishingly small coupling. [H, [36[ We extend this 
QMC method appropriately, allowing us to study spin- 
polarized superfluid neutron matter and therefore pro- 
vide the first benchmark calculation of this system using 
an ab initio microscopic simulation approach. 

In such systems, it is conventional to use the following 
measure of the population imbalance: 



P 



N^-N 



(1) 



where Nf and JVj, are the numbers of spin-up and spin- 
down particles, respectively, and P is called the polar- 
ization. The regime of large polarization is related to a 
question that has a long history in the framework of the 
BCS theory. In the BCS approach, superfluidity arises 
from the pairing of particles of different spin occupying 
states of opposite momenta near the Fermi surface. In the 
case of spin imbalance, the Fermi surfaces of the two com- 
ponents no longer coincide, making pairs with zero total 
momentum difficult to form. At some finite polarization, 
the gap between the two Fermi surfaces becomes so large 
that the system undergoes a quantum phase transition 
to a normal state (this is known as the Chandrasekhar- 
Clogston limit). 

Our aim in this work is to provide quantitatively re- 
liable results for superfluid low-density neutron matter. 
We therefore limit our simulations to small polarizations, 
P < 0.1, exploring the regime where pairing is likely to 
be energetically stable. We assume there are no exotic 
superfluid phases and no phase separation. We calculate 
ground-state energies at different total number densities 
(p = (N^+N±)/L 3 ), more specifically at p\ = 6.65 xl0~ 4 , 
p 2 = 2.16 x 1CT 3 , and p 3 = 5.32 x 1CT 3 fin" 3 . To put 
these densities into perspective we can compare them to 
nuclear matter saturation density: they are 0.41, 1.35, 
and 3.32 percent, respectively, of po = 0.16 fm -3 . At 
each total density, we study the cases of 35 + 33, 37 + 33, 
and 39 + 33 particles (see below) . We also compute the 
momentum distributions and pair-distribution functions 
for the two different components. 



II. QUANTUM MONTE CARLO 
A. Hamiltonian 

As pointed out in the Introduction, we do not need 
to include NNN interactions, since we are interested in a 
density regime where these are quite small. Thus, we use 
the following non-rclativistic Hamiltonian: 



H 



N 

£ 

k=l 



2m 



(2) 



where N = + N± is the total number of particles. 
The full neutron-neutron interaction is complicated, hav- 
ing onc-pion exchange at large distances, an intermediate 
range spin-dependent attraction by two-pion exchange, 
and a short-range repulsion. As already discussed, how- 
ever, in dilute neutron matter the dominant contributions 
come from the opposite-spin pairs, and specifically from 
the scattering length and the effective range, along with 
a short-range repulsive core which is important so as to 
avoid a collapse to a higher-density state 

In this work we are including an excess of neutrons of 
one species. Thus, it is also significant to take into ac- 
count the same-spin interactions in Eq. ([2]), which we do 
by using the interaction introduced in Ref. [34j . This 
interaction includes a propagator (see next subsection) 
in which all opposite-spin pairs interact through the 1 So 
channel of the Argonne vl8 (AV18) [37| potential, which 
fits s-wavc nuclcon-nuclcon scattering very well at both 
low- and high-energies. Thus, in what follows, for the 
purposes of the evolution the spins arc considered to be 
"frozen" , with the majority species being called f and the 
minority species being called -J,. We explicitly include the 
p-wave interactions in the same-spin pairs, and pertur- 
batively correct the S = l,Ms = pairs to the correct 
p-wave interaction. We use the Argonne vA 1 JAVA') po- 
tential to determine the p-wave interactions. [38j. 

Since we're studying neutrons, the AV4' interaction 
can be written as follows: 

V4,(r) = v c (r) +v rT {r)o- 1 ■ cr 2 , (3) 

which in the case of the 5 = (singlet) pairs gives: 

v s (r) = v c (r) - 3v a (r) . (4) 

In turn, the contribution from the 5* = 1 (triplet) pairs 
has the form: 



v P (r) = v c (r) + v a (r) 



(5) 



The same-spin potential contribution is small, but as the 
population imbalance increases the relative weights also 
change accordingly (see section lTlI A[) . While still keeping 
the potential of Eq. ^ in the propagator of our QMC 
method for the opposite-spin pairs, we have introduced 
a perturbative correction by writing Eq. ([3]) in terms of 
the Majorana exchange operator, which exchanges the 
positions leaving the spins unaffected: 

Vi(r) = v c (r) + v a (r)(-2P M - 1) (6) 
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B. Variational and Green's Function Monte Carlo 

The first step in our microscopic simulation is a Vari- 
ational Monte Carlo (VMC) calculation. Variational 
Monte Carlo is a relatively simple combination of clas- 
sical Monte Carlo and the variational (Rayleigh-Ritz) 
method; it was first used by McMillan in the 1960s. It is 
based on a variational trial wave function that is a 
reasonably good approximation of the true ground-state 
wave function. It contains variational parameters that 
in principle should allow one to approach the true wave 
function (see the next subsection for more details). A 
VMC calculation uses Monte Carlo integration to mini- 
mize the expectation value of the Hamiltonian: 



{H)VMC = /dR|*MR)|2 ^ Ea 



(7) 



thus optimizing the variational wave function tyy. The 
fact that it is relatively easy to perform a VMC simula- 
tion allows us to examine various possibilities in placing 
the excess particles in different momentum states. 

It is customary to use the output configurations of a 
Variational Monte Carlo calculation as input to a more 
extensive calculation using the method known as Green's 
Function Monte Carlo (GFMC). This method works by 
projecting out the exact, lowest-energy eigenstate <3/o 
from a trial (variational) wave function Vfy by treating 
the Schrodinger equation as a diffusion equation in imag- 
inary time r and stochastically evolving the variational 
wave function for a "sufficiently" long time. 

The evolution operator e~ lHt becomes e~ Hr in imag- 
inary time, commonly written as e ,—( H — E T)r ^ -^gj-g the 
Et is called the trial energy. Applying this operator to 
the variational wave function and expanding in terms of 
the complete set of eigenstates gives: 

*(t) = e-^-^^v = J2 a ^ iEt ~ ET)T ^* 

i 

= a Q e-( E °- ET)Ti S>o, limr^oo. (8) 

The GFMC technique is implemented by discretizing r 
and expressing the imaginary-time propagator as 



d -(H~E t )t 



J e -(H-E T )A T > (9) 



where r = nAr. If we now define the short-time Green's 
function by: 



G(R,R') = <R| e - (H - BT)Ar |R') 



where R is the configuration vector R = (ri, r2 . . . r/v) 
of 37V dimensions, then we can use it to calculate the 
evolved ^(r) starting from a set of VMC configurations. 
The short-time Green's function can be conveniently ap- 
proximated using the Trotter-Suzuki formula: 

G(R,R') « e- v ^^(R\e- TA ^\n')e- v ^">^e ETAT 

= e -(y(R) + y(R')-2£ T )Al (R|e _ TA r| R / ) (n) 



which is accurate to order (At) 2 . 

To avoid the fermion-sign problem, we impose what 
is known as the "fixed-node approximation" . A fixed- 
node simulation leads to a wave function that is the 
lowest-energy state with the same nodes as the trial wave 
function "fy. The resulting energy E is an upper bound 
to the true ground-state energy. Thus, if one chooses the 
variational wave function so that it includes a number 
of parameters, [HJ these parameters can be optimized 
to give the best approximation to the ground-state wave 
function (see next subsection). 



C. Trial wave function 

In these VMC and GFMC calculations there is a need 
to express the wave function of the system in terms of 
specific coordinate-space states. To this effect, we use 
a finite number N of particles with Born-von Karman 
(periodic) boundary conditions in a cubic box of volume 
L 3 , and N is chosen to be large enough so that the sys- 
tem can be assumed to be in the thermodynamic limit. 
For neutron matter, this was shown to be approximately 
66 particles in Ref. (Hf. Using a Cartesian coordinate 
system, the quantized plane waves e lk "' r will have mo- 
mentum vectors of the following discrete form: 



(12) 



where the n x ,n y ,n z are integers. The shell number / is 
defined such that / = n 2 + n 2 + n 2 . Thus, there is only 
1 possible combination of the n x ,n y ,n z that gives I = 
0, 6 combinations that produce / = 1, 12 combinations 
that lead to / = 2 and so on. Neutrons are spin one- 
half fcrmions, therefore, for equal populations for the two 
components the system has a closed-shell structure when 
N = 2,14,38,54,66,.... 

The simplest possible approximation (which, strictly 
speaking, is applicable only to the case of closed shells) 
that can be used for the input variational wave function 
is to describe the particles as being in a free Fermi gas. 
This approach assumes no correlations in the wave func- 
tion and is equivalent to having a product of two Slater 
determinants, one for spin-up and one for spin-down: 



$s(R) = D t D. 



(13) 



The single-particle states in the Slater determinants are 



(10) <t> n {r k ) 



Another choice for $(R), one which can also describe 
pairing, is the well-known BCS wave function $scs(R-) 
in its form for fixed particle number (which reduces to 
the Slater case under specific conditions). This choice is 
agreeable for both physical reasons (it reflects the fact 
that fermions with an attractive interaction can form 
Cooper pairs in the ground state) and mathematical rea- 
sons (unlike the Slater wave function, it has nodal sur- 
faces which can be varied so as to minimize the fixed-node 
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GFMC energy). Furthermore, a computationally appeal- 
ing aspect of this wave function is the fact that it can be 
written down as a determinant. 1391 



In this formalism, a general wave function with n pairs, 
u spin-up and d spin-down unpaired particles can be writ- 
ten as: 



J 



^bcs(R-) = A{[(f){r iv )...(t)(r nn ,)] [V^t(r„+i)...VM-( r «+«)] [4'H(r {n+1 y)...^ di (r {n+d y)] } 



(14) 



The unpaired particles are placed in ip^ and ipjl single-particle states. We can write this wave function as the 
determinant of an M x M matrix where M = n + u + d. In this work, we are interested in the case of a "gaplcss 
supcrfluid" which for polarized neutrons translates to 33 opposite-spin pairs along with an excess of unpaired particles 
of one species. When we have 2 extra spin-up particles, the corresponding matrix is written as follows: 



0(>"2,1') 0(7*2,2') 
V 0(7*35,1') 0(7*35,2') 



(7*1,33') -01f(ri) V>2t( r l) \ 
(7*2,33') "01t(l"2) ^2t( r 2) 



0(7-35,33' 
I 



(15) 



^lt( r 35) -02t( r 35) / 



The pairing function 4>{r) is a sum over the momenta 
compatible with the periodic boundary conditions. In 
the BCS theory the pairing function is: 



0(7-) 



(16) 



and here it is parametrized with a short- and long-range 
part as in Ref. [Hj]: 



(r) = P{r) + 



aie 

n, I<I C 



(17) 



We choose the single-particle states, i/j^, to be plane 
waves so as to ensure momentum conservation. We pick 
their momentum by checking values near the minimum 
(at each density) of the quasiparticle dispersion. The 
latter is calculated using the odd-even energy staggering: 



A = E(N + 1) 



1 



[E(N) + E(N + 2)] 



(18) 



where N is an even number of particles. At each den- 
sity, the minimum of the dispersion lies at a different 
momentum, pjij As already mentioned, we used VMC to 
place the particles at different momentum states. For 
very small polarizations, the minimum system energy is 
expected to be identical to the minimum of the disper- 
sion, which follows from adding only one extra particle. 
This is indeed the result we find, the only exceptions ap- 
pearing at density ps and particle numbers of 39 + 33 and 
higher (see below). 

In practice, we also include Jastrow (correlation) terms 
in the variational wave function: 

i^j i'¥=j' si' 

where the unprimed (primed) indices refer to spin-up 
(spin-down) particles. The Jastrow parts are taken from 



a lowest-order-constrained- variational method |40| calcu- 
lation described by a Schrodinger-like equation: 

-— V 2 f(r)+v(r)f(r) = A/(r) (20) 
m 

for the opposite-spin /(r) and 



-V 2 f P (r)+v(r)f P (r) 



2hr 



f P (r) = Xf P (r) (21) 



for the same-spin fp(r). Since the f(r) and fp(r) we 
use are nodeless, they do not affect the final result apart 
from reducing the statistical error. Since we are using 
the fixed-node approximation, we know that the result 
we obtain for one set of pairing function parameters in 
Eq. (|17[) will be an upper bound to the true ground-state 
energy of the system. The parameters are optimized in 
the full GFMC calculation as in previous works [HI, [33[ , 
providing the best possible nodal surface, in the sense of 
lowest fixed-node energy, for that form of trial function. 
As mentioned in section TlIBl this upper-bound property 
allows us to get as close as possible to the true ground- 
state energy of the spin-polarized supcrfluid system. 



III. RESULTS 

A. Equation of state 

We first address the energy of spin-polarized low- 
density neutron matter versus polarization. We have 
studied three total densities p\ = 6.65 x 10 -4 , pi = 
2.16 x 1CT 3 , and p 3 = 5.32 x 1(T 3 fm~ 3 . A smaller den- 
sity would correspond to neutron matter that is closer to 
the neutron star surface, which in turn implies a smaller 
magnetic field, and is thus less likely to be polarized. 



TABLE I: Results for the ground-state energy divided with 
the total number of particles for two species of neutrons in- 
teracting via the AV4' potential at p3 = 5.32 x 10" 3 fm -3 . 



iVt + ATj. 


EfM [McV] E n [keV] E u [keV] 


33+33 


2.133(1) 


11.4(1) 


11.3(1) 


35+33 


2.178(2) 


12.2(1) 


10.7(1) 


37+33 


2.230(2) 


12.7(1) 


10.0(1) 


39+33 


2.286(2) 


13.8(2) 


9.6(1) 



Also, lower density is more difficult to propagate in imag- 
inary time to a satisfactory accuracy level, given that 
lower density leads to larger inverse energy and there- 
fore longer propagation. Reversely, we do not study even 
larger densities because then we would have to include 
3-body interactions in our approach. For pure neutron 
matter this would imply, first, going away from firm ex- 
perimentally constrained interactions (three-neutron in- 
teractions are commonly fit to N = Z nuclei) and, sec- 
ond, the necessity of using spin-isospin dependent wave 
functions, therefore disallowing the use of approximately 
70 particles and thus the simulation of the thermody- 
namic limit in the framework of a variational ab initio ap- 
proach. Furthermore, larger densities would imply that 
the afore-mentioned perturbative correction in the prop- 
agator (see section III A[) would break down. For these 
calculations to be quantitatively reliable, the number of 
opposite-spin pairs should not stray too much from the 
"canonical" 33 + 33 case. Thus, for each of the three den- 
sities, we address the cases of 35 + 33, 37+33, and 39 + 33 
particles. For the case of the largest density we have also 
examined 41 + 33 and 43 + 33 particles (i.e. up to nearly 
double the polarization), finding the same overall trend. 

In most works on polarized neutron matter, the de- 
pendence of energy on polarization is taken to be 
quadratic. (23T4291 ] However, these works address dense 
neutron matter, in which the pairing has already 
reached gap closure, in the presence of a strong mag- 
netic field. In other words, most works in the literature 
study normal matter, in which an isolated spin-flip only 
impacts particles near their respective Fermi surfaces. In 
this work we study low-density polarized neutron mat- 
ter, implying that superfluidity plays an important role: 
flipping a spin is equivalent to breaking a pair. Thus, the 
pairing gap plays a decisive role in producing the polar- 
ized state. Our results for the energy versus polarization 
exhibit a linear trend: 



(22) 



similarly to the case of ultracold atomic gases at unitarity 
[HJ, where Quantum Monte Carlo calculations found a 
linear dependence of the energy on polarization at small 
population imbalances. 

In this connection, it is relevant to examine the in- 
teraction between (polarized) quasiparticles. Summing 
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FIG. 1: (color online) Ground-state energy per particle (in 
units of the free Fermi gas energy) for spin-polarized neutron 
matter. Shown are QMC results at three different total den- 
sities. The overall trend is linear, the slope depending on the 
density. 



single-particle excitation energies (taken from Ref. |34j |) 
for particles that are placed at the appropriate minimum 
momentum in Eq. (|12[) we find that the full Quantum 
Monte Carlo results are thus approximately reproduced. 
This implies that the quasiparticles are weakly interact- 
ing. 

The a(p) coefficients we have extracted from these 
results for the three densities are 0.37(5), 1.01(5), and 
1.84(9) McV, respectively. In Fig. [T]we show the energy 
per particle at different densities versus polarization. To 
facilitate comparison between results at different densi- 
ties we have divided the energy per particle with the en- 
ergy of a free Fermi gas at the same total density: 



Efg 



10 m 



(23) 



We notice that, just like in the case of unpolarized neu- 



tron matter 



34j . when the density increases the energy 
in units of Efg drops, but the rate of the drop is also de- 
creasing. The slight deviation from linear behavior at p\ 
stems from the afore-mentioned necessity to propagate 
up to longer imaginary times (at least by a factor of 3) 
in comparison to the other cases. This is also the reason 
why the results for the bigger systems at that density 
have larger error bars. 

Further results for the energy of polarized neutron 
matter are given in Table |U which refers to the largest 
density we have studied, p% = 5.32 x 10 -3 fm -3 . At this 
density we find the maximum value of the perturbative 
correction, which is 9 percent of the total energy. As 
should be expected, for the 33 + 33 system the energies 
of the ft an d the 44- interactions are identical (within 
statistical error). As we increase the polarization, there 
is a clear trend toward the increase of the tt relative im- 
portance, which for the 43 + 33 system becomes nearly 
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FIG. 2: (color online) Momentum distribution versus {k/kp) 2 
for the two different spins, for the case of 37 t an d 33 4- par- 
ticles (total density p-j, = 5.32 x 10 -3 fm~ 3 ), shown as circles 
and squares respectively, kp is taken here to refer to the 

total density, kp = (Sir 2 pzj 3 . The behavior exhibited is 
commonly referred to as "Fermi surface mismatch". 



FIG. 3: (color online) Same-spin pair-distribution function 
as a function of the distance divided with a measure of the 
average interparticle spacing, for the case of 37 t and 33 4 par- 
ticles (total density p3 = 5.32 x 10 -3 fm -3 ), shown as circles 
and squares respectively. Also given is the same-spin pair- 
distribution function for a non- interacting Fermi gas (line). 



double the corresponding 44- contribution. In all cases, 
the same-spin contribution to the energy is small, not 
growing to more than 1 percent of the total energy. 



B. Distribution functions 

We have also used GFMC to calculate distribution 
functions at p 3 = 5.32 x 1CT 3 fm" 3 for 37 t and 33 4 
particles: in contradistinction to the case of the energy, 
the results for these functions are not upper bounds to 
the true ground-state results, but they are expected to be 
accurate (the error being of second order in j^o — ^v\)- 
Starting with the momentum distribution, we first dis- 
cuss our expectations using mean-field BCS theory as a 
gui de. BCS is not quantitatively accurate in this regime, 
34J but can provide qualitative understanding. In BCS, 
the momentum distribution is given by the following ex- 
pression: 



n(k) = \ 



1 



E(k) 



(24) 



where £(k) = e(k) — pi, the chemical potential is pi and 
e(k) = is the single-particle energy of a particle with 
momentum k. The elementary quasi-particle excitations 
of the system have energy: 



E(k) = V£(k) 2 +A(k) 2 



(25) 



Overall, this is close to a step function for small gaps, but 
it changes considerably in the strong coupling regime. In 
general, the spread of the momentum distribution around 



p is approximately 2 A. At this density, the system ex- 
hibits a gap of A = 1.05(11) which as a fraction of the 
Fermi energy is A/Ep = 0.17(2) implying that there is 
no clearly defined Fermi surface. [H| Even so, in the case 
of spin-polarized Fermi gases it is customary to use the 
language of weak coupling and speak of a "Fermi surface 
mismatch" . This follows from the fact that the Fermi en- 
ergy is proportional to P^m and in this case the densities 
for the two spin populations are different. 

In Fig. [5] we show the momentum distribution com- 
puted using GFMC. This is calculated as the Fourier 
transform of the one-body density matrix, through: 



~ L 3 



U * v (ri,...,r n )J 



(26) 

where the curly brackets denote a stochastic integration 
over the angles. The integral over 5r = \r' n — r n \ is per- 
formed on a line analytically to avoid statistical errors 
due to the oscillatory radial dependence. In both cases, 
we see a considerable spread around the chemical poten- 
tial value, but we also notice a clear distinction in how 
the two species behave around that point, the majority 
species showing a "lag" in its decline. 

We have also computed the pair-distribution functions 
and have plotted them in Fig. [3] These are calculated 
from an expectation value of the form: 



gp(r)=AY,{*o\6(r ij -r)0? j \* v ) 



(27) 



%<3 



where we are interested in the case in which the oper- 
ator is simply unity, and the normalization factor A is 
such that gi(r) = g c (r) goes to one at large distances. 
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Such pair-distribution functions provide sum rules re- 
lated to density- and other response functions versus 
density and momentum. The solid line in the figure 
shows the pair-distribution function of noninteracting 
(NI) fermions with parallel spins: 

9c I (r) = l- [s\n{k F r)-k F rcos{k F r)f . (28) 

The x-axis is the interparticle distance divided with a 
quantity, r*o, which describes the average interparticle 
spacing: 

\*rl = - . (29) 
P 

The free Fermi gas result is close to but distinct from both 
QMC results (for spin-up and spin-down particles) due 
to the effect of the interactions. As is to be expected, the 
majority component values are slightly larger than those 
of the minority species, implying that it is slightly more 
likely to find a spin-up particle than a spin-down one. 

IV. CONCLUSIONS AND FUTURE WORK 

In summary, we have studied supcrfluid spin-polarized 
low-density neutron matter at small polarizations using a 
variationally optimized approach that includes the domi- 
nant well-known terms in the Hamiltonian. We have cal- 
culated the equation of state with the AV4' interaction at 
different densities. We have also calculated the momen- 
tum and pair-distribution functions. We find clear signals 
of a Fermi surface mismatch, as expected, and also a lin- 
ear dependence of the energy on the polarization. These 
results are in principle relevant to the physics of magne- 
tars. Furthermore, they could be tested directly by using 
ultracold fermionic atom gases with unequal spin popula- 
tions. In the case of cold atoms, Quantum Monte Carlo 
simulations of spin-polarized matter have been used as 



input to density- functional theory approaches. [4JJ, |42j 
Thus, our corresponding results for neutron matter might 
also be used as input to self-consistent mean-field models 
of nuclei. 

This line of Quantum Monte Carlo calculations, hav- 
ing first been applied to and verified in cold atomic ex- 
periments, can also provide directions for future work 
in the field of nuclconic infinite matter. The simplest 
case is that of a two-c omp onent gas where the two pop- 
ulations are equal. [H, [34j The next step is to examine 
the ramifications of taking different populations for the 
two components: this is the case of spin-polarized low- 
density neutrons studied in this work. Cold-atom ex- 
periments have by now also addressed Efimov physics, in 
which three components are involved. In the nuclear con- 
text, adding a third species could provide further insight 
into the physics of neutron stars. If the third component 
particles were taken to be protons and, as in this paper, 
only a few of them were added, then it would be possi- 
ble to study highly asymmetric nuclear matter. Another 
possible avenue of future research is related to optical lat- 
tice experiments with cold atoms: to first approximation 
these are equivalent to periodic external potentials. In 
the nuclear case, an external potential would allow us to 
study the static response of neutron matter and would 
also facilitate the understanding of the impact on neu- 
tron pairing of the ion lattice that exists in a neutron 
star crust. 
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